On population extinction risk in the aftermath of a catastrophic event 
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We investigate how a catastrophic event (modeled as a temporary fall of the reproduction rate) in- 
creases the extinction probability of an isolated self-regulated stochastic population. Using a variant 
of the Verhulst logistic model as an example, we combine the probability generating function tech- 
nique with an eikonal approximation to evaluate the exponentially large increase in the extinction 
probability caused by the catastrophe. This quantity is given by the eikonal action computed over 
"the optimal path" (instanton) of an effective classical Hamiltonian system with a time-dependent 
Hamiltonian. For a general catastrophe the eikonal equations can be solved numerically. For simple 
models of catastrophic events analytic solutions can be obtained. One such solution becomes quite 
simple close to the bifurcation point of the Verhulst model. The eikonal results for the increase in 
the extinction probability caused by a catastrophe agree well with numerical solutions of the master 
equation. 
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I. INTRODUCTION 

Evaluation of the extinction risk of a population in the 
aftermath of a catastrophe - a drastic deterioration of en- 
vironmental conditions - is of great importance in pop- 
ulation biology A practical approach to modeling of 
this problem assumes that (i) the population dynamics 
is governed by a (Markov) stochastic birth-death process 
0,13, (ii) this process has an absorbing state at zero pop- 
ulation size, and (iii) the catastrophe introduces an ex- 
plicit (and possibly strong) time dependence into one or 
more of the transition rates. The question we want to ad- 
dress in this work is the following: What is the extinction 
probability of the population at the time the catastrophic 
event is over and the environmental conditions return 
to normal? The explicit time dependence of the transi- 
tion rates, brought about by the catastrophe, makes the 
problem difficult for analysis. We will present here a for- 
malism that helps develop an insight into this class of 
problems. As a prototypical example of the population 
dynamics we will adopt a variant of the stochastic Ver- 
hulst logistic model. We will be interested in the regime 
of parameters where, if there is no catastrophe, a long- 
lived quasi-stationary distribution of the population sizes 
is formed. Because of the presence of an absorbing state 
at zero the population ultimately goes extinct in this 
model even without a catastrophe. A natural measure 
of the extinction dynamics is Vo(t): the probability that 
the population goes extinct by the time t. The focus of 
attention in this work is on the increase in the extinction 
probability Vo(t) caused by the catastrophe, see Fig. Q] 
We will assume that, if there is no catastrophe, the mean 
time to extinction (MTE), although finite, is too long 
to be relevant. Let Vq^) be the extinction probability 
before the catastrophe occurred. Neglecting a transient 
related to the relaxation to the quasi-stationary distribu- 
tion (and the exponentially small extinction probability 
during this stage), this probability grows very slowly ac- 



FIG. 1: (Color online) A schematic plot of the time-dependent 
extinction probability Vo{t) showing the effect of a catastro- 
phe with a characteristic duration T. 



is the MTE if there is no catastrophe. During the catas- 
trophe, which starts at t = t c , the extinction probabil- 
ity Vo(t) grows more rapidly and, after the catastrophe, 
reaches a value Vq . From then on (again, neglecting a re- 
laxation transient), the extinction probability grows very 
slowly, approximately as V^(t) ~ 1 — (1 — V^)e ( ~ tc ~ t ">/ r . 
We assume that the effect of the catastrophe is not too 
weak, "Pq ^> Vq , where Vq is the (exponentially small) 
accumulated extinction probability before the catastro- 
phe. The goal of this work is to evaluate the increase of 
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cording to Vq (f) ~ 1 - e' 1 ! 1 ', see e.g. 0, H 0. Here r 



the extinction probability AT^o = 

The formalism that we will employ to achieve this 
goal begins from a transformation of the master equa- 
tion (with time-dependent transition rates) into an exact 
evolution equation for the probability generating func- 
tion. For the problem in question this equation is a 
second-order linear partial differential equation (PDE). 
We then apply a time-dependent eikonal approximation 
to this evolution equation. This approximation, in many 
ways similar to the time-dependent WKB approximation 
in quantum theory Q, brings the second-order PDE to 
a first-order PDE of Hamilton-Jacobi type, with an ef- 
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fective classical Hamiltonian which explicitly depends on 
time. The further analysis deals with the characteristic 
lines of the Hamilton-Jacobi equations: the phase trajec- 
tories generated by the effective classical Hamiltonian. 

We will assume throughout the paper that both the 
initial population size, and the expected population size 
immediately after the catastrophic event, are much larger 
than unity. Under this assumption (and not too close to 
the bifurcation point of the Verhulst model, see below) 
we will evaluate the (exponentially large) net increase 
AVq in the extinction probability caused by the catas- 
trophe. This formalism yields an asymptotically correct 
value of the corresponding large exponent which makes 
it advantageous compared to the widely used van Kam- 
pen's system size expansion of the master equation and 
related methods, see e.g. 0, ||. Indeed, as it is well es- 
tablished by now, the Fokker-Planck equation, resulting 
from the van Kampen's system size expansion, cannot 
predict the MTE of stochastic birth-death systems with 
exponential accuracy. This can be traced to the failure of 
the Fokker-Planck equation in describing the actualprob- 
ability distribution tails of these systems 0, H, @, U Q ■ 

Here is the layout of the remainder of the paper. We 
begin Sec. [TT]with a brief reminder on the stochastic Ver- 
hulst logistic model. Then we derive the evolution equa- 
tion for the probability generating function and employ 
an eikonal approximation in order to determine the MTE 
of the population if there is no catastrophe. Here our re- 
sults coincide, with exponential accuracy, with those of 
previous works on the stochastic Verhulst logistic model. 
Section IHII deals with the effect of a catastrophic event 
on the population survival. We first demonstrate the ef- 
ficiency of the eikonal method by finding numerically the 
most probable path to extinction, and computing the cor- 
responding increase in the extinction probability due to 
the catastrophe, for a typical example of a catastrophic 
event. Then we obtain non-perturbative (in the catastro- 
phe magnitude) analytical results by adopting a simple 
schematic form of the catastrophic event: we postulate 
that the reproduction rate of the population drops in- 
stantaneously to zero at a specified time and recovers to 
the pre-catastrophe value, again instantaneously, after a 
given time T has elapsed. Section IIVI compares our nu- 
merical and analytical eikonal results for the extinction 
probability with direct numerical solutions of the master 
equation. A brief summary and discussion of our results 
is presented in Sec. |Vl 



II. VERHULST MODEL, PROBABILITY 
GENERATING FUNCTION AND EIKONAL 
APPROXIMATION 

As a prototypical example of self-regulating dynamics 
of an isolated population we consider a variant of the 
Verhulst logistic model: a Markov single-step stochas- 
tic birth-death process with continuous time. If there is 
no catastrophe, the reproduction and mortality rates are 



given by 



A„ = B n and fi r , 



Bn 2 



(1) 



respectively, where we use the same notations as in Ref. 
[9(. The reproduction rate per individual is constant, 
while the mortality rate per individual is constant at 
small population sizes, but grows proportionally to the 
population size when the latter is sufficiently large, ac- 
counting, for example, for competition for resources. For 
brevity, time and the transition rates in Eq. |T]) are 
rescaled with respect to the value of the mortality rate 
at small population sizes. 

At the level of deterministic modeling, the dynamics 
of the population size is described by the rate equation 



h(t) = (B — l)n(t) 



B 

N 



n 2 (t). 



(2) 



For B > 1 this equation has an attracting fixed point 
n s = N(l — 1/B) and a repelling fixed point n = 0. 
Throughout the paper we assume n s 3> 1; this necessarily 
requires N 3> 1. A linear stability analysis of Eq. $2$i 
around n — n s yields the characteristic relaxation time 
To = (B — 1) _1 (in the units of the mortality rate at small 
population sizes). 

Demographic stochasticity in the Verhulst model is ac- 
counted for by the master equation 



dV„ 

dt 



B{n - l)V n -i - BnPn 



B(n + 1) 
N 



V. 



n+l 



Bn 2 



where V n (t) is the probability that the population size at 
time t is n. The stochastic Verhulst model, as described 
by Eq. ([3]), and related models were considered in many 
works, see [i| [HI HH, EH E3, 13 an d references therein. It 
was found that, under certain conditions specified below, 
a long-lived quasi-stationary distribution, conditioned on 
non-extinction, is formed in this system after the relax- 
ation time 0(tq). The quasi-stationary distribution has 
a peak with relative width proportional to iV -1 / 2 around 
the attracting state n s of the deterministic description. 
At much longer times, however, the population goes ex- 
tinct. This is because n = is an absorbing state, so 
a rare sequence of events brings the process there with 
certainty. We will work in such a parameter regime that, 
if there is no catastrophe, the MTE is too long to be 
relevant. That is, we will be interested in times which, 
although much longer than the relaxation time To, are 
still much shorter than the MTE of the system without 
a catastrophe, that is the one described by Eq. ©. 
Notice that the same master equation ([3]) also describes 

the stochastic dynamics of three chemical reactions: A — > 
2A, A A 0, and 2A A [H|, with rates A = B , fx = 
1 + B/N, and a = 2(/i — 1), respectively. For definiteness 
we will use the notation of the Verhulst model in the 
following. 
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We will first demonstrate the method on the (well 
known) case of population extinction when there is no 
catastrophe. Introduce the probability generating func- 
tion @, Hi 



G(p,t) = J2p n V n (t), 



(4) 



where p is an auxiliary variable. G(p, t) encodes all the 
probabilities V n (t), as those are given by the coefficients 
of its Taylor expansion around p = 0. The generating 
function obeys the normalization condition 



G(l,t) = l 



(5) 



which follows from the conservation of the total probabil- 
ity. The distribution moments can be expressed through 
the p-dcrivatives of the generating function at p = 1, e.g. 
(n)(t)=j: n nV n (t) = d p G(p,t)\ p=v 

By multiplying Eq. ([3]) by p n and summing over all 
n one obtains, after some algebra, an evolution equation 
forG(p,t): 



dG B . . d 2 G 



R 1 B \ dG 



(6) 

The signature of the (empty) absorbing state is the 
absence of a term proportional to G(p,t) in Eq. ©. 
This immediately brings about the stationary solution 
G(p, t) = 1 corresponding to the empty state. 

In contrast to the Fokker- Planck equation, which is 
derivable from the master equation ([3]) by the van Kam- 
pen's system size expansion, the starting point of our 
theory - the evolution equation ((6j) - is exact. Through- 
out the rest of the paper we assume N 3> B. As our 
main results below are obtained with exponential accu- 
racy, this strong inequality enables us to neglect the term 
(B/N)dG/dp, and Eq. © becomes 



dG B . 8 2 G . 1WD . dG _ 

^(1-P)P ¥ + (P-1)(BP-1)^. (7) 



at 



This second-order PDE can be interpreted as an 
imaginary-time Schrodinger equation dG/dt = HG, 
where 



q2 

P)Pq^ + (P-1)(BP- I) 



dp 



is the Hamiltonian operator. 

In the framework of the spectral theory [f| 0, , the 
solution of the initial value problem for Eq. ([7]) can be ob- 
tained via an expansion in the eigenfunctions of a Sturm- 
Liouvillc problem related to the non-Hcrmitian operator 
H. The boundary conditions, needed for this Sturm- 
Liouville problem, come from the demand of bounded- 
ness of the eigenfunctions at the singular points of Ti. In 
this way one obtains a zero eigenvalue Eq = which de- 
scribes the true (absorbing) steady state of the system: 



extinction of the population. In addition, one obtains 
a discrete set of negative eigenvalues {E n }^ =1 that de- 
scribe the decay with time of the discrete set of eigen- 
functions. The relative weight of each eigenfunction in 
the expansion is determined, via G{p,t = 0), by the ini- 
tial probability distribution V n (t = 0). For N ^> 1 (and 
not too close to the bifurcation point B — 1) the higher 
eigenvalues E2, E%, . . . scale as Tq 1 , whereas the funda- 
mental eigenvalue E\ is exponentially small. Therefore, 
there are two widely separated time scales in the problem. 
During the fast relaxation time 0(tq) the higher modes 
decay exponentially, and the probability distribution ap- 
proaches the quasi-stationary distribution (QSD) men- 
tioned above, which is derivable from the fundamental 
eigenfunction. The exponentially small decay rate of the 
fundamental mode (and, correspondingly, of the time- 
dependent metastable distribution that is proportional 
to the QSD) is equal to the inverse MTE of the system. 
We checked that, as in other similar systems [a, ||a [l6| , the 
spectral theory yields accurate results both for the com- 
plete probability distribution and for the MTE, including 
the pre-exponent, in the catastrophe-free regime. How- 
ever, the spectral theory cannot, in general, accommo- 
date time-dependent transition rates: the main subject of 
this work. Therefore, in what follows we will use instead 
the time-dependent eikonal approximation [J] . Although 
it only gives exponential accuracy, the eikonal approxi- 
mation is readily generalizable to multi-population prob- 
lems [l7j and, as we will demonstrate shortly, to time- 
dependent Hamiltonians. 

Now we use the eikonal ansatz G(p, t) = exp[— <S(p, t)\ 
in the evolution equation ([7]). As can be checked a pos- 
teriori, S can be written as S = Ns(p, t), where s(p, t) = 
0(1). Therefore, the term (B/N)(l - p)pd pp S = 0(1) 
is small compared to the other terms which scale as N. 
Neglecting it, we arrive at a first-order PDE for S(p,t) 



dS 

~dt 



B 

N 



(1 



(p-l)(i?p-l) 



as 

dp 



(8) 

that can be interpreted as a Hamilton- Jacobi equation 
in the momentum representation, where p is the mo- 
mentum. Introducing the canonically conjugate coordi- 
nate q = —dS/dp, we obtain a one-dimensional classi- 
cal Hamiltonian flow with the time-independent Hamil- 
tonian 



^ ( J^PI 



Bp+l)q 



(9) 



The fact that H(l,q) = reflects the probability con- 
servation, see Eq. ([5]), whereas the fact that H(p, 0) = 
manifests the presence of the absorbing empty state. 

It is convenient to shift the momentum p = p — 1, leav- 
ing the coordinate q unchanged. The new Hamiltonian 
is 



H(p, q)=P 



N 



l)g + B(p+l)-l 



9, (10) 
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FIG. 2: (Color online) The phase plane (p, q) emerging in 
the eikonal approximation to the stochastic Verhulst model. 
The fixed points are denoted by M (the mean-field point), 
F (the fluctuational point) and (the trivial point). Solid 
lines depict zero-energy trajectories, including the heteroclinic 
trajectory MF. Dashed lines show two examples of non-zero- 
energy trajectories. 



while the Hamilton equations are 

* = l^ = i2p+1) { Bq -§ q2 )- q > 

9H . 2 . / 2B \ 
P = -^ = -(p 2 +P)[B- W q)+p. (12) 

This Hamiltonian flow was investigated, in the context of 
the three above-mentioned chemical reactions, by Elgart 
and Kamenev [l5| . The phase plane (p, q) , defined by the 
Hamiltonian (I10p . provides a useful visualization of the 
extinction dynamics, see Fig. [2j As the Hamiltonian does 
not depend explicitly on time, it is an integral of motion: 
H(p, q) = E = const, and the problem is integrable. The 
zero energy trajectories, E = 0, play a special role here, 
as will become clear shortly. One type of zero-energy 
trajectories are the mean field trajectories, staying on 
the line p = 0. Here Eq. (JTTJ) becomes 



q=(B-l)q-^q 2 



(13) 



which coincides with the rate equation @. This fact 
provides the interpretation of the q- variable as the reac- 
tion coordinate. Importantly, the attracting fixed point 
of the mean-field, or deterministic line p = 0, q m f = n s — 
N(l — 1/B) becomes a hyperbolic point (0, q m f) of the 
phase plane (p, q) . We call this fixed point the mean- field 
point and denote it by M. 

The Hamiltonian (|10p has two more zero-energy lines. 
One of them is the extinction line q = which includes 
two more hyperbolic fixed points of the phase plane (p, q): 
the "fluctuational point" (1/B — 1, 0), denoted by F, and 
the trivial point (0,0) denoted by 0. The third zero- 
energy curve 



includes a heteroclinic trajectory that exits, at time t = 
— oo, the mean field point M along its unstable eigendi- 
rection, and enters, at time t = oo, the fluctuational 
point F along its stable eigendirection. This heteroclinic 
trajectory represents the optimal (most probable) path, 
or the instanton connection, of the extinction dynam- 
ics. Indeed, it describes, in the eikonal language, the 
most probable sequence of discrete events bringing the 
system from its quasi-stationary state to extinction, cf. 
[1 [H, [H, HI . If there is no catastrophe, the MTE of the 
population is, with exponential accuracy, r ~ exp(iSo), 
where the (zero-energy) action is 



So 



qpdt . 



(15) 



and the integration should be performed along the zero- 
energy heteroclinic trajectory ([33]). The result is 



r ~ exp 



Qo(p)dp 



exp N 



B-l-lnB 



B 



(16) 

The approximation is valid when So 3> 1. One can check 
that the result (1161) coincides, up to a pre-exponent, with 
those of [IdO, El S3. 

Of a special interest is the region of parameters close 
to the bifurcation point of the Verhulst model: TV -1 / 2 <c 
B — 1 -c 1, where the left inequality is required for the 
validity of the eikonal approximation. For B — 1 <C 1 the 
action 



So 



N(B-l) 



(17) 



scales as the square of the distance to the bifurcation 
point. It has been found recently [Tt], S3] that, close to 
the bifurcation point, the Verhulst model, the SIS (Sus- 
ceptible - Infected - Susceptible) model of epidemiology 
[1 Gill! OS OS 13, the SIS model with demogra- 
phy the SIR ( Sus cep tible - Infected -Recovered) 
model with demography [TtI SI SH and other related 
stochastic models become, in the leading order, identical 
if their rates are properly rescaled. One can check that, 
in this limit, the term —(B/N)pq in the square brackets 
in Eq. (fTU|) can be neglected compared to the rest of the 
terms, and one arrives at the "universal" Hamiltonian 



H (p, q)=P 



' N 



p + B-l) q 



(18) 



introduced in Ref. [l5( in the context of the three chem- 
ical reactions mentioned above. All three zero-energy 
lines in the phase plane of the universal Hamiltonian arc 
straight lines, and the action So, given by Eq. (fT7|) . is 
simply the area of the triangle formed by these straight 
lines. 
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III. CATASTROPHIC EVENT AND ACTION 
CALCULATION 

To model a catastrophic event we assume that, because 
of unfavorable environmental changes, the reproduction 
rate drops and then recovers to the pre-catastrophe value. 
This can be described by introducing a time-dependent 
factor f(t), such that /(±oo) = 1, into the reproduction 
rate: 



Xn(t) =£/(*) n : 



(19) 



The mortality rate ji n remains constant in time in this 
model, see Eq. (JlJ. At the level of deterministic model- 
ing, the population size is described by the rate equation 



h(t) = [Bf(t)-l] n{t)-^n 2 (t), 



(20) 



the master equation becomes 



dt 



Bf(t)(n-l)V n -i-Bf(t)nV„ 



1 



B(n + l) 
N 



V n +1 - n 



Bn* 



Vn, (21) 



while the evolution equation for G(p,t) is again an 
imaginary-time Schrodinger equation dG/dt = TLG with 
the Hamiltonian operator 



H 



#a-p)^+(p-i)P/(«)p-< 



that now explicitly depends on time. The same 
eikonal ansatz G(p,t) = exp[— <S(p, t)] brings about the 
Hamilton- Jacobi equation for S(p, t) which defines a clas- 
sical Hamiltonian flow with the time-dependent Hamil- 
tonian 



H(p,q,t) =p 



B 

N 



(p+l) 9 + B/(t)(p+l)-l 



q, 



with the same coordinate q and shifted momentum p = 
p — 1 as in the time-independent case. The Hamilton 
equations are 



q = ^ = (2P + D 



P 



dp 

dH , 2 x 



Bf(t)q~^q 2 



Bf{t) §q 



q, (23) 
+ p. (24) 



As the Hamiltonian H is not an integral of motion any- 
more, the problem is in general non-integrable. There 
are, however, two planes in the extended phase space 
(Pi q-i t) where the Hamiltonian is still conserved and equal 
to zero. These are the mean-field plane {p = 0, q, t) and 
the extinction plane (p, q = 0, t). In the mean- field plane 
Eq. (jTTJ) becomes 



q = [Bf(t) - 1] q 



B 



(25) 



which again coincides with the rate equation (|20[). In the 
eikonal approximation the extinction occurs along the 
instanton connection of the time-dependent Hamiltonian 
(f2"2"| in the extended phase space (p, q, t) . To reason- 
ably describe a catastrophe the function f(t) should be 
bounded from above. In this case there exists a hetero- 
clinic trajectory [p op (t), q op (t)] which exits the mean-field 
fixed point M well before the catastrophe and arrives at 
the fluctuational fixed point F after the catastrophe is 
over. The full action along this heteroclinic trajectory 



S = 



{-Qop(t)pop(t) - H[q op (t),p op (t),t}}dt (26) 



[where H is given by Eq. (|22|) ]. determines (the logarithm 
of) the increase in the extinction probability caused by 
the catastrophe. 



A. Instanton calculations: numerical solution 

For a smooth f(t), the Hamilton equations (f2"5)l and 
([24| (and similar Hamilton equations for other models) 
can be accurately solved by a shooting method with a sin- 
gle shooting parameter. As in the unperturbed case, the 
instanton connection must exit, at t = -co, the mean- 
field point M and enter, at t = oo, the fluctuational 
point F. In a numerical solution we specify, at some 
initial moment of time t = ti n prior to the catastrophe, 
the coordinate q(U n ) and momentum p(U n ) < which 
lie on the unperturbed heteroclinic trajectory (114|) in a 
vicinity of the mean-field point M. The initial momen- 
tum p{ti n ) can play the role of the shooting parameter. 
It must be chosen sufficiently close to zero for the nu- 
merically found instanton to approximate well the en- 
tire perturbed instanton. On the other hand, the instan- 
ton must reach, at a final time tf after the catastrophe, 
(a close vicinity of) the fluctuational point F. There- 
fore, for too a small p(ti n ), intrinsic logarithmic slow- 
down near the fixed point tt increases the computation 
time and causes accumulation of numerical errors. De- 
noting the numerically found instanton (or optimal path) 

by \php (t), 3op (t)]) we can write the action along it as 

Snum = {-q$(t)p%\t)-H [q$>(t),p%>(t),t]}dt. 

(27) 

As an example, let us consider 



/(*) 



AB (t-ta) 2 
ft 

B 



(28) 



where AB < B is the catastrophe magnitude as mani- 
fested in the change of the reproduction rate, t c is the 
time when the catastrophe reaches its maximum, and T 
is the catastrophe duration. Figure [3] shows projections 
on the (p, q) plane of the numerically found instantons 
for several values of the catastrophe duration T. Wc 
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FIG. 3: (Color online) Projections on the (p, q) plane of the 
"extinction instantons" found by solving the Hamilton equa- 
tions l(23)l and ((24)) numerically for fit) given by Eq. J25J. 
The parameters are N = 200, B = 2, AB = 0.75, and 
t c — tin = 40. The catastrophe durations are T = 1 (solid 
line 1), T = 3 (solid line 2) and T = 20 (solid line 3). The up- 
per dashed line is the instanton when there is no catastrophe 
and f(t) = 1 = const. The lower dashed line is the instanton 
for the system with B replaced by B — AB — const. The up- 
per and lower dashed lines yield the upper and lower bounds, 
respectively, on the action as a function of T. 



compared the numerically found action (|27|) to (the log- 
arithm of) the increase of the extinction probability due 
to the catastrophe, AVq, determined by solving numer- 
ically the master equation (j!H]) with fit) from Eq. 
Such a comparison is presented in Sec. 4. 



B. Instanton calculations: analytically soluble 
example 

We present now an example when the eikonal equa- 
tions, describing the impact of a catastrophe on the 
population extinction probability, are exactly soluble. 
We assume that the reproduction rate A„ drops instan- 
taneously to zero at t = and recovers to the pre- 
catastrophe value, again instantaneously, after time T 
has elapsed. In terms of the function fit) we have 



. 1 if t < or t > T , 
JW - i o if < t < T. 



(29) 



see Fig. |4j With this f{t) the solution of the determin- 
istic rate equation ([2!>]) is 

( [B(e* - lRlf 1 , 0<t<7\ 
nit) = n s x < 

I [B(e T - l) e -(B-W-T) + 1] - 1 , t > T , 

(30) 

where n(t = 0) = n s = N(l — 1/B) corresponds to the 
mean-field point M. The deterministic solution, shown 
in Fig. predicts a complete recovery of the population 
after the catastrophe. However, the stochastic effects, 




FIG. 4: (Color online) A model catastrophic event for which 
the eikonal problem is exactly soluble. The catastrophe starts 
at time t = 0, when the reproduction rate drops to zero, and 
ends at t = T, when the reproduction rate recovers to its 
pre-catastrophe value. 




FIG. 5: (Color online) Effect of the catastrophic event, de- 
scribed by Eq. (|29[) . on the expected population size as pre- 
dicted by the solution (|30|l of the deterministic rate equation 
(|25|) . The parameters are N — 200, -8 = 2, and T = 1. 



missed by the rate equation, can be greatly enhanced 
because of the temporary decline in the population size. 
As a result, they can increase the extinction probability 
considerably. As we rely on the eikonal approximation 
in describing this effect, we are interested in the regime 
where, in spite of the catastrophe, the expected popu- 
lation size by the end of the catastrophe, n(T), is still 
sufficiently large: n(T) 3> 1. 

Because of the special shape of the function f(t), there 
are now two distinct Hamiltonians: the unperturbed 
Hamiltonian (|10p before and after the catastrophe and 
the zero-reproduction-rate Hamiltonian during the catas- 
trophe: 



H c (p,q) 



-P 



(p+l)q + l 



(31) 



Each of the two Hamiltonians is an integral of motion on 
the corresponding time interval. The instanton can be 
found by matching three separate trajectory segments: 
the pre-catastrophe, catastrophe and post-catastrophe 
segments. Figure [5] shows a projection of the instanton 
on the (p, q) plane. To recall, the instanton must exit, at 
t = — oo, the mean-field point M and enter, at t = oo, the 



7 



fluctuational point F. The matching conditions at times 
t = and t = T are provided by the continuity of the 
functions q(t) and p(t). The pre- and post-catastrophe 
segments must have a zero energy, E — 0, so they are 
parts of the original zero-energy trajectory of the un- 
perturbed problem, see Eq. (fl4|) . For the catastrophe 
segment, however, the energy E = E c is non-zero and a 
priori unknown. It parametrizes the intersection points 
Pi and pi (see Fig. [6]) between the unperturbed zero- 
energy line 



qo (p) = N 



N 



B{l+p) ' 
and the non-zero-energy line H c = E c : 



Qc(p, E c 



N 



2B(l+p) 



4E C (1 



Np 



(32) 



(33) 



Solving the algebraic equation qo(p) = q c {p, E c ) for p, we 
obtain 



Pi{E c ) 



P2(E C ) 




(34) 



To determine E c we demand the duration of the catas- 
trophe to be T. Using Eq. for f(t) = and Eq. ([531). 
we obtain an algebraic equation for £J C = E C (T): 



dp 



MEc) p{(2B/N)(p+l)q c (p,E c ) + l 



(35) 



where q c (p,E c ) is given by Eq. (|3"3"|) . The net increase 
of the extinction probability of the population due to 
the catastrophe is proportional to exp[— S(T)], where the 
action S(T) is 



PllEa(T)] 



rP2[E c (T)] 
S(T) = / q (p)dp 
Jl/B-l 

+ I q {p)dp~E c (T)T 
Jpi[E c (T)] 

We can also rewrite this action as 
rPilEa(T)] 



q c [p,E c (T)]dp 



(36) 



S(T) = So- {q (p) - q c \p, E C {T))} dp 

Jp2[E c {T)\ 

- E C (T)T, (37) 

where Sq is the unperturbed action, see Eqs. (fT5|) and 
(|16p. and the integral term in this equation corresponds 
to the area of the shaded region in Fig. [Sj 

To obtain more visual results, let us consider the limit 
of B — 1 <C 1. In this limit the pre- and post-catastrophe 
Hamiltonian reduces to the universal Hamiltonian (|18[) . 




1/B-1 



FIG. 6: (Color online) Projection on the (p, q) plane of the 
"extinction instanton" (the thick solid line going from M to 
F) for the catastrophic event described by Eq. (129[) . Points p\ 
and p2 correspond to times t = and t = T where the catas- 
trophic event begins and ends, respectively. At t < and 
t > T the instanton follows the zero-energy heteroclinic tra- 
jectory q — qo(p), whereas at < t < T it follows a non-zero- 
energy trajectory q — q c (p,E c ). The energy E c — E C (T) is 
determined by the given catastrophe duration T. The area of 
the shaded region corresponds to the integral term in Eq. (|37|) . 



Furthermore, as the terms {B/N)pq and (B/N)q can 
be neglected here compared to 1 in Eq. (f3"Tj) . the zero- 
reproduction-rate Hamiltonian during the catastrophe 
simplifies drastically: 



H c (p, q) - -pq 



(38) 



and the catastrophe segment of the instanton trajectory 
becomes simply 



q c lp,E c ) 



E c 



(39) 



The corresponding Hamilton equation for p [Eq. I24p ]. on 
the catastrophe segment, is p ~ p, so ln(p2/pi) — T and 



^ ~ exp(T) 
Pi 



(40) 



In their turn, p\ — p\[E c ) and P2 — paC^e) are the roots 
of the equation qoljp) — q c lj>, E c ) which now reads 



NIB-l+p) ~ . 

P 



Equations ((40j) and (gTJ) yield 



1 - B 

Pi - T - , P2 

e 1 + 1 



(1 - B)e T 
e T + l ' 



and 



N(B-l) 2 , 2 (T 

E c ~ — i '- cosh" 2 - 

4 V 2 



(41) 



(42) 



(43) 
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One can see that, as T — > 0, the intersection points p\ 
and P2 merge signaling no change of the unperturbed 
action. As T — > oo, p\ approaches the mean-field point 
M, whereas P2 approaches the fluctuational point F. In 
the latter case the population quickly dies out. 

Now we use Eq. ([37jl to calculate the extinction action 
caused by the catastrophe. After a simple algebra we 
obtain 

S(T) ~ J* , (44) 
1 + e 1 

where the unperturbed action So is given by Eq. (|17[) . 
For short catastrophe durations, T < 1, we obtain a 
small correction, linear in T, to the unperturbed action: 

S(T)~S (l-T/2). (45) 

For long catastrophes, T> 1, the total action <S(T) de- 
cays exponentially with an increase of T, 

5(T) ~ 2S e~ T , (46) 

signaling a rapid extinction. 

Interestingly, one can express the extinction action 
(I44p in terms of an effective universal Hamiltonian that 
is catastrophe-free. Indeed, Eq. (|44|) can be interpreted 
as the area of a right-angled triangle: 

S(T)^^±n eff (T), (47) 

where {B — 1), the absolute value of p at the fluctuation 
point F (for B — 1 -C 1), is one of the legs of the triangle. 
The other leg, n e //(T), is the harmonic mean of n s and 
n(T): 



n e //(T) n s n(T) ' 

[To recall, close to the bifurcation n s ~ iV(B — 1) is the 
steady-state pre- and post-catastrophe population size 
and n(T) ~ n s exp(— T) < n s is the population size im- 
mediately after the catastrophe, as predicted by the de- 
terministic rate equation.] For short catastrophes n(T) 
is close to n s , and one obtains a small correction to the 
unperturbed action. For long catastrophes n(T) <C n Sl 
and n e ff(T) ~ 2n(T). We are unaware of a simpler way 
to arrive at these results via arguments based on the de- 
terministic description of the catastrophe, Eq. (|30[) . and 
on the knowledge of the unperturbed action, Eq. (Tl"T|) . 

For the eikonal approximation to be valid, we have to 
demand that both the total action, and the correction to 
it caused by the catastrophe, be much larger than unity. 
This yields a range of rescaled catastrophe durations 

5 " 1 <T<ln<S (49) 

which becomes broader as iSo 3> 1 increases. 

The eikonal theory provides, with exponential accu- 
racy [26j , the value of the increase of the extinction prob- 
ability AV = Vq - Vq (see Section HJ) in terms of the 
"catastrophe action" S(T): 

AV oc e- 5 ^ . (50) 




FIG. 7: (Color online) The extinction probability Vo(t) vs. 
time if there is no catastrophe. The parameters are N = 
10,800 and B = 1.08. Inset: ln[l - Po(t)] vs. time. The 
slope of this graph yields the MTE: r ~ 6.2 x 10 13 , so that 
lnr ~ 31.8. For comparison, the eikonal theory, see Eq. (|16|l . 
predicts lnr ~ 30.4. 



IV. COMPARISON TO NUMERICAL 
SOLUTIONS OF THE MASTER EQUATION 

We tested the predictions of the eikonal theory by solv- 
ing the (truncated) master equation ([3]) numerically. We 
observed that, if there is no catastrophe, the system ap- 
proaches a quasi-stationary probability distribution after 
a time 0(tq), as expected. This probability distribu- 
tion then very slowly decays in time, while the extinction 
probability Vo(t) very slowly grows. The corresponding 
decay/growth time is in good agreement with the MTE 
t predicted from Eq. (fl"6|) , see Fig. [7j 

We modeled the catastrophe by multiplying the unper- 
turbed reproduction rate by f(t), with f(t) from either 
Eq. pg]). or Eq. The validity of the eikonal ap- 

proximation demands So 3> 1. [One must also demand 
that the expected population size in the quasi-stationary 
state be much larger than unity: q m f = N(B — 1)/B ^> 1. 
However, above the bifurcation point B > 1, when the 
deterministic version of the Verhulst model has a non- 
trivial attracting fixed point, the criterion So 3> 1 is al- 
ways more restrictive.] One more criterion for the validity 
of the eikonal approximation is S(T) 3> 1, for all values 
of T that we used. 

In our numerical solutions of the master equation the 
initial condition corresponded to a fixed population size: 
V n (t = 0) = 5 n ,n , where <5n,n is the Kronecker delta, 
and no ^> 1. In this case an immediate extinction be- 
fore relaxing to the quasi-stationary state has an expo- 
nentially small probability. The catastrophe time t c was 
chosen to be several times longer than the relaxation time 
To for the chosen parameters, but much shorter than the 
expected MTE. Finally, the catastrophe duration T was 
chosen such that the correction to the action caused by 
the catastrophe and the total action satisfy the conditions 
So— S(T) ^> 1 and S(T) ^> 1. In the exactly soluble case, 
presented at the end of the previous section, these two 
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FIG. 8: (Color online) The extinction probability Vo(t) vs. 
time (at times much smaller than the MTE t), found by solv- 
ing numerically the (truncated) master equation (|2ip with 
f(t) from Eq. ([29]). (a): N = 14,400, B = 1.08 and T = 2.5. 
(b): N — 10, 800, the rest of parameters is the same as in (a). 
The starting time of the catastrophe in both cases, t c = 300, 
obeys the condition t r <C t c <C T. Vo(t) before the catastro- 
phe is negligibly small and cannot be seen in this scale. 




FIG. 9: (Color online) (a) A comparison between the natural 
logarithm of the net contribution AVo of the catastrophe to 
the extinction probability Vo(t) (see text), and —S(T) from 
Eq. (J2TJ) vs. T. (b) The ratio of the two quantities vs. T. 
The parameters are N = 200 and B = 2, AB = 0.75, and 

tc - tin = 40. 



inequalities reduce to the double inequality ([4"9"]1 . Figure 
[8] presents, for two different sets of parameters, our nu- 
merical results for Vo(t). Here f{t) is given by Eq. 

Figures 191 and [TUl compare the predictions of our eikonal 
theory with the numerical solutions of the master equa- 
tion. To estimate the direct impact of the catastrophe 
on the extinction probability we calculated the quantity 
APo = T-'o — ~Pq, where Vq is the measured extinction 
probability before the catastrophe starts but after the re- 
laxation stage ends. As t c <C r, most of the contribution 



FIG. 10: (Color online) The natural logarithm of the numeri- 
cally computed AVo vs. T (solid line) is compared to —S(T) 
predicted by Eq. (144 ti (dashed line). The parameters in (a) 
and (c) are the same as in Fig. [8] (a) and (b), respectively. 
In (b) and (d) the ratios of the same quantities, InAPo and 
—S(T), are plotted vs. T for the same parameters as in (a) 
and (c), respectively. The analytical result in (b) deviates 
from the numerical result by 2.8% at the maximum. The 
corresponding deviation in (d) is 4.4% indicating that, as So 
increases, the agreement between the analytical and numeri- 
cal results improves. 



to Vq comes, for the parameters and initial conditions 
we worked with, from the projections of the initial condi- 
tion on the ^rajjidly decaying) higher eigenmodes of the 
system, see 



Our numerical results for fit) from Eq. ([28]) are pre- 
sented in Fig. [HI The figure compares mAT^T 1 ) and 
—S(T) found by solving the eikonal equations numeri- 
cally, see Section 3.1. Figure [10] corresponds to the sim- 
ple catastrophe described by Eq. (|29p. Here In AVq(T) is 
compared to —S(T) for the two sets of parameters used 
in Fig. [8] These two sets of parameters were chosen to 
obey the strong inequality B — 1 <C 1, so that we could 
test the analytical prediction (|4"4"]) . Good agreement be- 
tween the theory and numerical computations is observed 
in all cases, over a broad range of T. One can see from 
Fig. [9] (b) and Figs. \W\ (b) and (d) that the agreement 
is best for intermediate values of T and deteriorates for 
small and large values of T . This behavior is consistent 
with the validity criteria of the eikonal approximation: 
for too a small T the inequality So — S(T) 3> 1 does not 
hold, whereas for too a large T the inequality S(T) ^> 1 
does not hold. The quantity — (In AVo) /S includes in- 
formation about the pre-exponent A(T) of the complete 
relation AVq — A(T)e~ s ( T \ This pre-exponent has not 
yet been determined theoretically. 
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FIG. 11: (Color online) Not fully reversible catastrophes, as 
reflected in the reproduction rate vs. time, and the extinction 
instantons for B\ > B2 (a and c) and B\ < B2 (b and d). 
The solid lines in (c) and (d) are the instantons, the dashed 
lines are the zero-energy lines, see Eq. (|14J1 . for B = B\ and 
B = B 2 . 



V. SUMMARY AND DISCUSSION 

We have employed a model example of the stochas- 
tic Verhulst logistic model to show that a catastrophic 
event (modeled as a temporary drop of the reproduc- 
tion rate) exponentially increases the extinction proba- 
bility of an isolated self-regulated stochastic population. 
This was achieved by combining the probability gener- 
ating function technique with an eikonal approximation, 
and by calculating the eikonal action over the instanton 
of the effective (time-dependent) classical Hamiltonian of 
the system. We have shown that eikonal approximation 
to the evolution equation for the probability generating 
function provides a robust and efficient way of evaluat- 
ing the impact of catastrophic changes, reflected in the 
reproduction/mortality rates, on the stochastic popula- 
tion dynamics. Once the corresponding large parame- 
ter is present in the problem, the eikonal approximation 
can be applied to a broad class of stochastic population 
model, single- or multiple-species. 

Although we dealt with fully reversible catastrophes, 
the eikonal method is easily extendable to more com- 
plicated situations. For instance, consider the case when 
the reproduction rate B reaches, after a catastrophe ends, 
a value B 2 different from the pre-catastrophe value B\. 
The extinction instantons, and the corresponding total 
actions, can be easily found for B\ > B 2 and for B\ < B 2 , 
see Fig. QTJ This immediately yields, with exponential 
accuracy, the corresponding increase in the extinction 
probability caused by the catastrophe. 

The theory presented in this paper is non-perturbative 



in the catastrophe magnitude. For relatively weak and/or 
short catastrophes one can develop an eikonal perturba- 
tion theory by assuming that the correction to the unper- 
turbed action So is small compared to the unperturbed 
action itself (but still much larger than unity, so that 
the eikonal theory is valid). For a time-periodic modula- 
tion of the reaction rates such a theory has been recently 
developed in Refs. [27|, [28|. We checked that, for the ex- 
actly soluble catastrophe model, presented at the end of 
section Hill such a perturbation theory correctly predicts 
the small- T asymptote of the action, Eq. (gSj) . 

Let us put the results of this work into a more gen- 
eral context of employing eikonal approximation for a de- 
scription of large fluctuations in stochastic populations, 
the reproduction and/or mortality rates of which depend 
on time, reflecting different types of environmental vari- 
ations. In the present paper we have dealt with the 
case of a catastrophe, where the reproduction rate under- 
goes a strong temporary drop. The papers [27], l28| em- 
ployed eikonal approximation to deal with rate modula- 
tions which are periodic in time (modeling daily, monthly, 
or annual cycles of environmental parameters) . The peri- 
odicity of the rate modulation makes it possible to apply 
special theoretical tools unavailable otherwise (such as 
the Kapitsa method, applicable when the rate modula- 
tion period is very short compared to the relaxation time 
towards the long-lived quasi-stationary state of the pop- 
ulation (28|). Finally, a variant of eikonal approximation 
has been recently applied to the case where the reproduc- 
tion and/or mortality rates are modulated stochastically: 
by noisy environmental variations [29j | . The stochastic 
environment may be thought of as a random sequence 
of weak and strong catastrophic events, with a certain 
probability distribution of their amplitudes and dura- 
tions. This probability distribution is sensitive to the 
color of the environmental noise. Understanding the im- 
pact of a single catastrophe, achieved in the present work, 
turns out to be vital in elucidating the role of the noise 
color in the population extinction[29(. 

The main conclusion of this work is that eikonal ap- 
proximation yields a simple, accurate and visual means 
of assessing the viability of stochastic populations which 
experience catastrophic changes of environmental condi- 
tions. 
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